9 随机区组试验设计

9.1 加载环境与数据准备

9.1.1 安装包检测并自动安装（必须联网）

pkgs <- c("tidyverse","agricolae","here","knitr","latex2exp","cowplot")                   
not_installed <- pkgs[!(pkgs %in% installed.packages()[ , "Package"])]   
if(length(not_installed)) install.packages(not_installed)  

9.1.2 加载R必备包

library(tidyverse)
library(agricolae)
library(here)
library(knitr)
library(latex2exp)
library(cowplot)
options(digits = 4)

9.1.3 试验数据读取

dt <- readxl::read_excel(here("RES/dataset/产量干物质.xlsx"),
                         sheet = "干物质",
                         skip = 2) 
dt <- dt %>% filter(Sampling == 1) %>% select(-Sampling)
dt

表 9.1 试验数据集

PlotID	Rep	Irrigation	Biochar	Grainyield	DMA
1	R1	CF	B0	8.578	18.75
2	R1	CF	B20	8.653	16.98
3	R1	CF	B20M	7.852	15.79
4	R1	AWD	B0	7.811	16.16
5	R1	AWD	B20	8.794	17.56
6	R1	AWD	B20M	6.509	14.05
7	R2	CF	B0	7.474	16.93
8	R2	CF	B20	7.806	14.21
9	R2	CF	B20M	7.386	15.71
10	R2	AWD	B0	6.956	14.58
11	R2	AWD	B20	6.953	13.01
12	R2	AWD	B20M	7.878	16.60
13	R3	CF	B0	6.346	12.45
14	R3	CF	B20	8.881	18.98
15	R3	CF	B20M	8.157	17.87
16	R3	AWD	B0	9.809	19.93
17	R3	AWD	B20	7.234	13.80
18	R3	AWD	B20M	6.190	11.49

9.2 试验数据检测

方差分析数据必须满足3个基本条件，即独立、正态，方差齐次性。

9.2.1 方差正态性检测

Shapiro-Wilk检验用来检验数据是否符合正态分布，类似于线性回归的方法一样，是检验其于回归曲线的残差。该方法作者推荐在样本量很小的时候使用，比如N<20。但是也有作者推荐在大数据集上使用。该作者将这种修改后的方法运用在R语言的stats包中的shapiro.test函数中。该检验原假设为H0：数据集符合正态分布。

dt[c("Grainyield","DMA")] %>%  map_df(
  function(x){
    result <- shapiro.test(x)
    return(list(W=result$statistic,PValue=result$p.value))
  },.id = "Respone"
) %>% select('指标'=Respone,'W统计量'=W,'P值'=PValue)

表 9.2 产量数据和干物质数据正态性检测

指标	W统计量	P值
Grainyield	0.9749	0.8825
DMA	0.9813	0.9625

由表 9.2中P值列可知，所有指标shapiro方差正态性检测P值均大于0.05，表明结果接受了原假设，所有指标均符合正态分布。

9.2.2 方差齐次性检验

R中最常用的方差齐性检验（bartlett检验），bartlett.test()调用格式为bartlett.test()调用格式为bartlett.test(x,g…)其中，参数x是数据向量或列表(list);g是因子向量，如果x是列表则忽略g，当使用数据集时，也通过formula调用函数。

双因素试验必须对两个因素均进行齐次性检测

9.2.2.1 对第1个因子（灌溉模式）进行齐次性检测

首先，检验第1因子（灌溉模式）不同水平下方差是否等值

dt[c("Grainyield", "DMA")] %>%  map_df(function(y) {
  trt <- dt$Irrigation
  result <- bartlett.test(y ~ trt)
  list(
    '统计量' = result$statistic,
    '自由度' = result$parameter,
    'P值' = result$p.value
  )
}, .id = "指标") %>%
  kable(digits=3)

表 9.3 产量数据和干物质数据第1个因子的齐次性检测

指标	统计量	自由度	P值
Grainyield	1.069	1	0.301
DMA	0.275	1	0.600

9.2.2.2 对第2个因素子（生物炭）进行齐次性检测

接着，检验第2因子（生物炭）不同水平下方差是否等值

dt[c("Grainyield", "DMA")] %>%  map_df(function(y) {
  trt <- dt$Biochar
  result <- bartlett.test(y ~ trt)
  list(
    '统计量' = result$statistic,
    '自由度' = result$parameter,
    'P值' = result$p.value
  )
}, .id = "指标") %>%
  kable(digits=3)

表 9.4 产量数据和干物质数据第2个因子的齐次性检测

指标	统计量	自由度	P值
Grainyield	1.067	2	0.587
DMA	0.197	2	0.906

由表 9.3和表 9.4可知，产量和干物质的两个因子bartlett.test检测P值均大于0.05，表明结果接受了原假设，所有指标在两个因子上均符合齐次性检测。

9.3 方差分析

由此可见，试验数据满足正态和方差齐次性，可正常进行方差分析。

以表 15.1数据集为例，采用随机区组试验设计模型对试验数据进行方差分析，结果如表 10.2。

response <- dt$Grainyield
ANOVA <- aov(response ~ Rep + Irrigation * Biochar + Error(Rep), data = dt) %>%
  summary()
ANOVA <- pmap_df(list(data=ANOVA,name=names(ANOVA)),
        ~.x[[1]]) 
ANOVA %>%
  kable(digits=3)

表 9.5 随机区组试验方差分析表

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Rep	2	1.177	0.589	NA	NA
Irrigation	1	0.499	0.499	0.493	0.498
Biochar	2	1.650	0.825	0.815	0.470
Irrigation:Biochar	2	2.542	1.271	1.256	0.326
Residuals	10	10.120	1.012	NA	NA

有的学者可能对方差贡献率感兴趣，利用各因子及交互效应等平方和数据，可进一步得到各因子和交互效应贡献率情况，如表 10.3。

ANOVA %>% mutate(Contri=`Sum Sq`/sum(`Sum Sq`)*100) %>%
  kable(digits=3)

表 9.6 方差分析+方差贡献率

	Df	Sum Sq	Mean Sq	F value	Pr(>F)	Contri
Rep	2	1.177	0.589	NA	NA	7.363
Irrigation	1	0.499	0.499	0.493	0.498	3.123
Biochar	2	1.650	0.825	0.815	0.470	10.321
Irrigation:Biochar	2	2.542	1.271	1.256	0.326	15.899
Residuals	10	10.120	1.012	NA	NA	63.294

9.4 事后均值对比

9.4.1 定义HSD结果提取函数

getHSD <- function(hsd){
  merge(x=hsd$means %>% select(std),
        y=hsd$groups,
        by="row.names") %>% 
    rename(trt=Row.names)
}

9.4.2 对第1因子进行HSD均值检测

DFerror <- ANOVA["Residuals         ","Df"]
MSerror <- ANOVA["Residuals         ","Mean Sq"]

means <- dt$Grainyield
hsd.Irrigation <- HSD.test(
  means,
  trt = dt$Irrigation,
  DFerror = DFerror,
  MSerror = MSerror
) %>%
  getHSD() %>% select(trt, means, std, groups)

hsd.Irrigation %>%  kable(digits=3)

表 9.7 不同灌溉模式下产量均值对比结果

trt	means	std	groups
AWD	7.570	1.149	a
CF	7.904	0.785	a

9.4.3 对第2因子进行HSD均值检测

hsd.Biochar <- HSD.test(means,
                        trt = dt$Biochar,
                        DFerror = DFerror,
                        MSerror = MSerror) %>%
  getHSD() %>% select(trt, means, std, groups)

hsd.Biochar %>%  kable(digits=3)

表 9.8 不同生物炭水平下产量均值对比结果

trt	means	std	groups
B0	7.829	1.230	a
B20	8.053	0.841	a
B20M	7.329	0.804	a

9.4.4 对交互因子进行HSD均值检测

hsd.Interaction <- HSD.test(
  means,
  trt = interaction(dt$Irrigation, dt$Biochar),
  DFerror = DFerror,
  MSerror = MSerror
) %>%
  getHSD() %>% select(trt, means, std, groups) %>% separate(col = trt, into =
                                                              c("Irrigation", "Biochar"))

hsd.Interaction%>%  kable(digits=3)

表 9.9 不同灌溉模式和生物炭水平下产量均值对比结果

Irrigation	Biochar	means	std	groups
AWD	B0	8.192	1.464	a
AWD	B20	7.660	0.992	a
AWD	B20M	6.859	0.897	a
CF	B0	7.466	1.116	a
CF	B20	8.446	0.566	a
CF	B20M	7.799	0.388	a

9.5 出图

9.5.1 设置试验因子各水平排序

hsd.Irrigation$trt <- hsd.Irrigation$trt %>% fct_relevel("CF", "AWD")
hsd.Biochar$trt <- hsd.Biochar$trt %>% fct_relevel("B0", "B20", "B20M")
hsd.Interaction$Irrigation <-  hsd.Interaction$Irrigation %>% fct_relevel("CF", "AWD")
hsd.Interaction$Biochar <- hsd.Interaction$Biochar %>% fct_relevel("B0", "B20", "B20M")

mylabels.irr <- c(CF = TeX("$I_{CF}$"), AWD = TeX("$I_{AWD}$"))
mylabels.bio <- c(
  B0 = TeX("$B_{0}$"),
  B20 = TeX("$B_{20}$"),
  B20M = TeX("$B_{20M}$")
)

9.5.2 第1因子灌溉模式不同水平均值对比

G.irr <- ggplot(hsd.Irrigation ,aes(x=trt,y=means,fill=trt ))+
  geom_bar(stat="identity")+
  geom_errorbar(aes(ymin=means-std,ymax=means+std),
                width=0.2)+
  geom_text(aes(y=means+std,label=groups),vjust=-0.2)+
  cowplot::theme_cowplot(font_size = 8,line_size = 0.4)+
  scale_fill_discrete(labels=mylabels.irr )+
  scale_x_discrete(labels=mylabels.irr )+
  ylim(0,12)+
  labs(x = "Irrigation Regime",
       y = TeX("Grain yield ($t~ha^{-1}$)"),
       fill = "Irrigation")
G.irr

图 9.1 不同灌溉模式下产量情况

9.5.3 第2因子生物炭不同水平均值对比

G.Bio <- ggplot(hsd.Biochar ,aes(x=trt,y=means,fill=trt ))+
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = means - std, ymax = means + std), width = 0.2) +
  geom_text(aes(y = means + std, label = groups), vjust = -0.2) +
  cowplot::theme_cowplot(font_size = 8, line_size = 0.4) +
  scale_fill_discrete(labels = mylabels.bio) +
  scale_x_discrete(labels = mylabels.bio) +
  ylim(0, 12) +
  labs(
    x = TeX("Biochar application rate ($t~ha^{-1}$)"),
    y = TeX("Grain yield ($t~ha^{-1}$)"),
    fill = "Biochar"
  ) 
G.Bio

图 9.2 不同生物炭水平下产量情况

9.5.4 交互因子不同处理间平均值对比

dod <- position_dodge(width = 0.95)
G.int <- ggplot(hsd.Interaction , aes(x = Irrigation, y = means, fill = Biochar)) +
  geom_bar(stat = "identity", position = dod) +
  geom_errorbar(aes(ymin = means - std, ymax = means + std),
                position = dod,
                width = 0.2) +
  geom_text(aes(y = means + std, label = groups),
            vjust = -0.2,
            position = dod) +
  cowplot::theme_cowplot(font_size = 8, line_size = 0.4) +
  scale_x_discrete(labels = mylabels.irr) +
  scale_fill_discrete(labels = mylabels.bio) +
  ylim(0, 12) +
  labs(
    x = TeX("Biochar application rate ($t~ha^{-1}$)"),
    y = TeX("Grain yield ($t~ha^{-1}$)")
  )
G.int

图 9.3 不同灌溉模式和生物炭水平下产量情况

9.5.5 合成并出图

GG <- plot_grid(
  plot_grid(G.irr, G.Bio, nrow = 1, labels = letters),
  G.int,
  nrow = 2,
  labels = c("", "c")
)
GG

图 9.4 主效应和交互效应